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ABSTRACT 

Context. Compton scattering is involved in many astrophysical situations. It is well known and has been studied in detail 
for the past fifty years. Exact formulae for the different cross sections are often complex, and essentially asymptotic 
expressions have been used in the past. Numerical capabilities have now developed to a point where they enable the 
direct use of exact formulae in sophisticated codes that deal with all kinds of interactions in plasmas. Although the 
numerical computation of the Compton cross section is simple in principle, its practical evaluation is often prone to 
accuracy issues. These can be severe in some astrophysical situations but are often not addressed properly. 
Aims. In this paper we investigate numerical issues related to the computation of the Compton scattering contribution 
to the time evolution of interacting photon and particle populations. 

Methods. An exact form of the isotropic Compton cross section free of numerical cancellations is derived. Its accuracy 
is investigated and compared to other formulae. Then, several methods to solve the kinetic equations using this cross 
section are studied. 

Results. The regimes where existing cross sections can be evaluated numerically are given. We find that the cross section 
derived here allows for accurate and fast numerical evaluation for any photon and electron energy. The most efficient 
way to solve the kinetic equations is a method combining a direct integration of the cross section over the photon and 
particle distributions and a Fokker-Planck approximation. Expressions describing this combination are given. 
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Introduction 

. Compton scattering has become a major process in astro- 
' physics, and for more than half a century, it has been ap- 
plied to a large variety of astrophysical situations, particu- 
larly in high energy, low density media. It was first intro- 

■ duced in the context of e lectron cosmic rays interacting with 
, thermal stellar photo ns (jFollinlll947t iFeenberg fc Primakofil 

■ 119481: lDonah^ll95lL etc.). It was then used to model gases 
' of high energy electrons in the strong radiatio n field of X- 
, ray or relativistic sources (jLevich fc Svunvaevlll97lh . It is 

■ also responsible for the up-scattering of low energy photons 
' from the microwave background radiation ([Zeldovich et al.l 

Il972t l. More recently, much effort has been spent on high 
energy media in relativistic sources such as 7-ray bursts, 
AGN, and X-ray binaries. Particular focus is now on the 
synchrotron self-Compton radiation that is thought to play 
a major role for exa mple in AGN ([Tavecchio et al.ll200ll : 
ISauge fc He'^l2004D . 

The many varied applications result in different energy 
regimes for the photons and the particles. The energy of 
particles spans from the very cold electrons of pulsars winds 
to the hot gas in the corona of X-ray binaries (Ew 100 keV, 
i.e. E/irieC^ ~ 1), and to high energy cosmic rays (up to 
E > 100 TeV, i.e. E/rrieC^ > 10*). The energy of photons 
also spans a wide range: in the case of AGN for instance, 
it goes from low energy synchrotron photons (A = Im, i.e. 
hv/rrieC^ « 10~^^) to 7-ray emission (with E > 10 TeV, i.e. 
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hu/meC^ > 10^ for some sources). General formulae must 
be able to deal simultaneously with very small and very 
large energies. 

The exact cross section is quite complex and does not 
give any physical insight into the fundamental process. 
Moreover, often, only a fraction of the particle and pho- 
ton spectra contribute significantly to the Compton scat- 
tering in one given astrophysical situation. For these rea- 
sons most of the theoretical work has been done on de- 
riving limiting and/or averaged forms that can be easily 
included in the modeling of astrophysical objects. For ex- 
ample, much work is now based on the diffusion approach 
for particles, valid when large angle scattering can be ne- 
glected, that is, when the energy of the particles does not 
va ry significantly in one single interaction (work initiated 
bv lKompaneetsl[l957f ) . Analytical expressions averaged over 
Maxwellian distributions and thus only valid in purely ther- 
mal plasmas have also been used. Many other approxi- 
mations have been proposed in general papers, r eviews 
and books (see ^umenthal fc Gould 1970; Zcldovic hjUfTSt 
'Rv bickifcLighmanlll979l : lKat^ll987l : INagirner fc PoutanenI 
1993, for example). Nevertheless the relevant energy can- 
not always be confined to a narrow range where simple 
approximations can be made. Also more general expres- 
sions that can be used regardless of a specific problem 
are required for generic codes designed to address vari- 
ous situations. In these cases an exact description with- 
out limitation on the photon and electrons energies is re- 
quired. A few exact formulae have been proposed for the 
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isotropic Compton cross section ([jonesi Il968t iBrinkmamil 
11984 iNagirner fc PoutanenI 11994 NP94 hereafter). Such 
formulae are analytically equivalent. When evaluated nu- 
merically however, they behave differently with resp ect to 
accura cy issues. For example, the original formula by I Jonesi 
(|l968l ) fails to describe the up-scattering of low energy 
photons which represents a typical astrophysical situation. 
Although these problems are well known, the Jones' cross 
section is still widely used in the literature. 

Such general expressions are not very helpful for physi- 
cal understanding. Nevertheless, they are complementary 
to asymptotic forms for they provides quantitative re- 
sults that can be compared to observations. It has be- 
come particularly true recently, since more and more 
codes have been developed to deal with photon-photon, 
photon-particle and p article-partic l e interactions in high 
energ y plasmas (e.g. Coppi 119921 : Poutanen fc SyenssonI 
19961: [Navakshin fc Melial Il998l: Belmont et al.l l2008t 



Vurm fc PoutanenI 12009^ " Even with an accurate, exact 



cross section, numerical computation of the Compton con- 
tribution to the evolution of the particle and photon dis- 
tributions is problematic. When the distributions and cross 
section are discretized in energy bins, the grid resolution 
puts severe constraints on the code accuracy. Typically, 
high energy particles are scattered down by numerous 
and inefficient scattering events off soft photons. If the 
resolution is insufficient, these individual events are not 
captured by the numerical scheme and the global evo- 
lution of the particle dist r ibutio n is not described accu- 
rately. iNavakshin fc MeTial (|1998^ proposed that a combi- 
nation of this method with a Fokker-Planck approxima- 
tion could enable good acc uracy. This idea was applied in 
a recent numerical code by lBelmont et al.l (|2008l ) (see also 
I Vurm fc Poutanenll2009| ). Here we investigate the method 
and quantify its accuracy. 

In the first section, we derive a new form of the distri- 
bution of scattered photons. The accuracy of this formula 
is discussed and compared to other expressions. In the sec- 
ond section, numerical errors of several kinetic methods are 
estimated and compared. 

1. Evaluating the isotropic distribution of 
Compton-scattered photons 

1.1. Exact spectrum 

As long as stimulated scattering is not considered, the evo- 
lution of a photon population interacting by Compton scat- 
tering with an electr on distribution is describ ed by the fol- 
lowing equation (see lRvbickifcLighmajilll979l for example) : 



dpofeiPo) J dko X 



dt 

U{kQ)c^{po,ko ^ k) - /^(fc)c-^(po,fe -> feo) 
ak dko 



(1) 



Here — dMuj/ (Pr / (fik and /e — dAfe/d'^r/d^p are 
the photon and lepton distributions in their respective 6- 
dimension momentum space. The first term in the brackets 
describes the contribution of all photons of momentum fco 
being scattered to momentum k after one interaction. The 
second one describes photons of momentum k being scat- 
tered to any other momentum. Both are proportional to 
the Compton differential cross section da/dk{po, ko fc) 



giving the probability of one photon of momentum fco being 
scattered to a momentum fc by an electron of momentum 
Pq. This cross section can be interpreted as the distribu- 
tion of scattered photons resulting from the interaction of 
mono-energetic photons with one single electron and it will 
be often be referred to as such hereafter. 

When the photon and particle distributions are 
isotropic, this equation can be integrated over all direc- 
tions. By noting d^p/(mec)^ = p^dpdflp and d^ k / {mec)^ = 
Lu'^dujdfti.j (where oj = hv IrrieC^), it yields: 



dt 



da , 
dio 



uj)dNjLOo)dNM 



-N^{lu) J cao{po,ujo)dNe{po) (2) 

where N^{uj) = Attui'^ f^{k) and Np{p) — 47rp^/p(p) are the 
angle-integrated distributions. Here 



is the total scattering cross section, and 

dCle.o f dQ,i^ of da 



da 



An 



An 



dwdd,. 



-dn^ 



(3) 



(4) 



is the angle-averaged cross section for isotropic Compton 
scattering. 




Fig. 1. Angles and notations. The incoming photon direc- 
tion no is defined by the incident angle ^o and an az- 
imuthal angle ^po that orientates the global picture (not 
shown here). The scattered photon direction n is defined 
by the scattering angle and the azimuthal angle (j). In 
coordinates (61,62,63), f3o/Po — (Mo,sin0OiO) and n — 
(cq , sin a cos (j), sin 0) , so that fi — iioCa + sin a sin fio cos (j). 



Eq.Ujcan actually be simplified further. Because of sym- 
metry, the differential cross section da / du / dVl^ does not 
depend independently on the direction of the incoming pho- 
ton and lepton, but only on the angle between their direc- 
tions: /xo = cos6'o, 6*0 = (/3o,no) where (3q = Po/jo is the 
lepton velocity in units of the speed of light and no is the 
direction of the incoming photon (see Fig.[T]). Then, Eq. |3| 
reduces to: 



da . 



d^o 
2 



da 



dujdfl. 



-dn^ 



(5) 
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Also, although the result of a single scattering event is de- 
scribed by 6 variables (energy and direction of the scattered 
photon and particle), these variables are not independent. 
The conservation of the total energy-momentum 4-vector 
sets 4 constrains, so that only 2 independent variables are 
necessary to describe the scattering outcome. Although 
other choices can be made, these are often chosen to be 
the angles describing the scattered photon direction n: the 
cosine of the scattering angle Cq. = cos a, a = {nt),n) and 
an azimuthal angle (see Fig. [T]). Then, the differential 
cross section da / duj / dflc^ is redundant and includes a delta- 
function that ensures that the total energy-momentum is 
conserved: 



da 



.,4 



da 



dcn 



-5{uj — uj)d^od4>dco 



Here the latter condition is given by: 



UJ = UJQ 



1 - /3oMo 



l-/3oA* + ^(l-c, 



(6) 



(7) 



where fj, — /iqCq + — /i§ — cos is the cosine of the 
angle between the direction of the scattered photon and the 
incoming electron. 

The angle differential cross section is a well known result 
of quantum mechanics. In the rest frame of the electron, the 
general Klein-Nishina cros s section does not depend on the 
azim uthal angle and reads (|Heitleij|l954l : iRvbickifcLighmajil 
fl979D : 



■'dc' 



1 da' 
2^dZ 



1 

1 

P' 



(8) 



where ax is the Thomson cross section, p'(cq) — i^'/^o ^ 
1 / {1 + ujq{1 — c'^)) is the ratio of the scattered to incoming 
photon energy. Here, all variables denoted with ' refer to 
quantities in the rest frame of the particle. 

Eg. [6l ari d m can be shown to be equivalent to Eq. 12 of 
iJonesI (fl968l) and Eq. 6.2.1, 8.1.1 and 8.1.3 of NP94. Getting 
the angle-averaged cross section is now quite straightfor- 
ward. However, including the Klein-Nishina cross section 
in Eq. [S] requires changes of frame and the triple integral 
quickly leads to cumbersome expressions. Different integra- 
tion variables can be used and the triple integral can be per- 
for med in sever al orders. Here we follow the original method 
bv lJonesI (|l968t) . All quantities are expressed in the electron 
rest frame which gives: 



da 
du! 



da' S[u! — uj] d/iQ d(j)' 
d<:7o'(l + /3oMo)'~2^ 



(9) 



Integration is then performed over fj.'^, and then cf)'. 
The latter integration is made with the vari able 1 -I- 
/3oA*'(</'')- The resuhing formula bv I JonesI (|l968l ) contains 
several misprintfQ. Using slightly different notations, the 
correct formula for the scattered distribution reads (see also 
iPe'er fc Waxmanll2005D : 



da 
du! 



3(Tt 



870P0W0 



(10) 



^ Contrary to what is claimed in ICoppi fc BlandfordI (|l990l ). 
all intermediate expressions are correct and only the final equa- 
tions 23-27 contains misprints. 



with 



F 




(11) 



4- {j/uj + "fo/ujQ) 

- — 4 + (70^0 -I- jus) 

X- y 

acosh ^ ■\/ A+a;+/2^ / \/A+ ^ 

asin (^|A_|x_/2^ /VTM if A_ < 
asinh (^|A_|x_/2j / ^/\>Z\ if A_ > 

where 70 and 7 = 70 -t- ti^o ~ ^'^^ the Lorentz factors 
of the incoming and scattered particles respectively and 
P = Y^7^ — 1 is the momentum of the scattered particle. 
Here the following notations have been used: 



A+ 
A_ 

c± 



= pI + 270^^0 + ^0 

= - 270a; + O)^ 



= 1 + UJLJ() + 2\± 



(12) 
(13) 
(14) 
(15) 



The cross section is proportional to the difference [F]^ of 
function F evaluated at the upper and lower boundaries of 
the last integration (& and a respectively) . At these bound- 
aries, the variables x± have the following expressions: 



min [(70 + Pq)/ujo; (7 +p)/^^] 
min [(70 +po)/^ ; (7 +P)/^o] 



1/{u;luox1116) 

l/iiULOox'LJU) 



The critical photon frequencies for which the arguments of 
the min and max functions are equal define distinct do- 
mains of the scattered distribution. For the lower integra- 
tion boundary a, both arguments are equal at a unique fre- 
quency: Lo = ujQ. For the higher boundary, both arguments 
are equal at two critical frequencies: u; = ujq and uj — ujc, 
where 



^0 



70 +Po 



70 - Po + 2^0 



(18) 



The scattered distribution is obviously continuous at these 
frequencies. However, higher derivatives are not, as it can 
be seen in Fig. 3] for example. 

The scattered distribution is bounded in energy and the 
maximal and m inimal energ ies guarantee that all radicals 
remain real (see [jMleslfT96a for more details): 



Wo 



70 -Po 



where 



1 



UJr, 



I0+P0 + 2a;o 

Wo + 70 - 1 if ^^0 > t^o(Po) 
LOc if loq < t^o(Po) 



(1 +Po - 70) 



(19) 
(20) 

(21) 



When, ujQ > ujq, the photon can gain up to the entire elec- 
tron energy. However, this absolute limit is not necessarily 
reached. When loq < ujq (as for Compton up-scattering for 
instance), only a fraction of the electron energy is at most 
transfered to the photon. 
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1.2. Cancelation issues 



As was first discussed by I Jone j ()1968[ ). evaluation of Eq. 
1121 suffers from accuracy issues due to machine round-off 
errors. Since most astrophysical problems involve an energy 
range spanning many orders of magnitude, Eg. 1121 includes 
combinations of very small and very large terms. In this 
form, most large terms must cancel out, which obviously 
represents a numerical challenge. An example of accuracy 
issue is shown in Fig. 3] 

Cancelation errors appear when the relative difference 
of two terms is smaller than the machine precision. Here we 
focus on computations in double precision, i.e. when reals 
are coded in 8 octets (64 bits) and their mantissa is coded 
in 52 bits. The relative error due to machine precision is 
then of the order of e = 2"^^ k. 10"^^. 

We have checked the accuracy of the scattered distri- 
bution for a large range of photon and lepton energies 
and identify the domain of the (wqiPo) plane where Eq. 
[T^ gives accurate results. Results are shown in Fig. [21 We 




Fig. 2. Accuracy domains. Equation [T^] is accurate for pho- 
tons energy and electron momentum (wqiPo) in the region 
delimited by th e solid line. The origina l formula derived 
by [J ones' flOGSf) , and the cross section by iPe'er fc Waxmai] 
(|200 5). have basically the same properties and are accu- 
rately evaluated in the same region. The published for- 
mula by NP94 is accurate in the domain determined by 
the dashed line. These boundaries are for double precision 
computing. The star shows the point for which the scat- 
tered distribution is plotted in Fig. [H 



find that Compton scattering of mild photons (wq 1) by 
mid-relativistic particles (po ~ 1) is well described by this 
formula but problems arise for lower or higher energies. 
In particular, the scattered distribution cannot be evalu- 
ated in the most typical astrophysical case, namely the up- 
scattering of soft photons by high energy particles. 

In the latter case, [Jones! ([l968) derived asymptotic ex- 
pansions in luq/jo << 1. However, a large number of orders 
must be kept when lvq /70 ^ 1 and the numerical evaluation 
often becomes time consuming. Moreover, other situations 
are also affected by numerical issues, for which such an 
expansion is not relevant. Here, we rather concentrate on 
analytical manipulation of the original formula to get an 
exact expression free of cancelation issues. 



1.3. Re-writing the exact cross section 

Numerical accuracy issues result from two kinds of cancela- 
tions: when the relative difference of function F evaluated 
at the two integration boundaries is very small and when 
some terms constituting F should cancel out. Here we deal 
with both successively. 



1.3.1. Cancelation issues when \F'' - F""] « F^.F^ 

This happens typically when for example, at 

the distribution upper and lower photon energies Wmin and 
^max where the distribution vanishes, but where F"" and F^ 
can remain large. 

Differences of hyperbolic functions between the upper 
and lower boundaries can be computed analytically by us- 
ing the following trigonometric relations: 

acosh(v^) — acosh(\/6) = asinh (^\J a{b — 1) — ■\/6(a — 1)^ 

asinh(v^) — asinh (Vfo) = asinh ^■\/ a(l + 5) — ^5(1 + a)j 

asin (\/a) — asin (Vfo) = asin (^y/ a{l — 6) — \/b{l — a)j 

The difference of the two hyperbolic functions between 
the two integration boundaries a and b simply reads: 



asinh(A-\/A±) if A± > 0, or asin(A-^— A±) if A± < 0, where 



A = L+(4)-L+«) 

= min [p, Po] + min [0; oj + ujq — p — po] /2 . 



(22) 

(23) 
(24) 



Although the last term of Eq. [T2| is not divergent it is 
numerically ill-behaved when A_ << 1 and it is best to use 
the following auxiliary function: 



Six) 



1 / asinh(^y^) 



asm{-\/—x) 




(25) 



which is evaluated as a series expansion for small argument 
(2n-(- l)!(-x)" 



Six) 



E 



^4"(n!)2(2?i + 3) 



for \x\ < 1 



(26) 



Differences of all other terms in Eg. [T2l between the two 
integration boundaries can also be computed analytically to 
factorize A and get the following quite simple expression for 
the scattered distribution: 



da _ 3A 

duj 8uj1^oPo 

with 



G 



(27) 



G 



(1 -I- uujo) 



z+ + z_ 



2 

LOLOq 



2[A]+- A2(1 + cjcjo) [zA\\ 



z^'^S (AA2 



(28) 



where A± = (l/z± + AiA^)-!/^ and z± = x%x\. The scat- 
tered distribution is proportional to A. This factor holds as 
the difference between the two integration boundaries a and 
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b and vanishes at the minimal and maximal photon energies 
ijJmin and Wniax- For high particle energy, it is best computed 
as: 

lo+l + Po+P 



A = min<min(p;po 



X min(a;; o^o) — 



2(p + Po) 

max(cij; luq) 



(29) 



Here we note that an accurate evaluation of the scat- 
tered electron momentum p is crucial. When the latter 
is low (typically when p < 10~^), the simple derivation 
p — \J (70 + lliq — a-')^ — 1 is not accurate enough and an 
alternative derivation must used: 



P^Po 



(cj -u;o)(p + Po)(7 + 7o) 
2po(p + Po) - (w - wo)(7 + 70) 



(30) 



Equations l??!!^ are accurate over a larger domain than 
the original formula and can be used in many astrophysical 
situations. 



1.3.2. Cancelation when 



z+ — z_| << z+, z_ 



However, this new expression reveals other differences hid- 
den in the original one and that result in cancelation is- 
sues, namely when \z+ — Z-\ << Z-. In such cases, the 
differences must be computed analytically. Although alter- 
nate expressions for the differences can often be accurate 
in more cases, it is sufficient to use them when the relative 
differences are smaller than e = 10^^°. 

By using the definition of and z_ , the differences in 
the first three terms of Eq. [25]can be computed as: 



2 

LULOq 

1 1 

UJ UJQ 



1, 



(7 -I- p -I- 70 +po)'^Lou;o 
ilo + Po) {7 + p){p + Pay 



[^A]t = '-±±^ [A]t + ^±±^ [zt 



respectively, where 
[A]+ = (cj + cjo)(7 + 7o) 



z r = - 



1 1 



(31) 

(32) 
(33) 

(34) 
(35) 



(7-1-^-1-70 +Po)^ 



1 1 



(p + Po)(7+P)(7o +Po)' w Wo 
In the same manner, the last term can be computed as: 



(36) 



cz'/'S 
where 



3/2 , 3/2 

Z I ' + z_ 



^3/2 



I 3/2 , 3/2 

+ C+Zi +c^z_ 



c+ + c- + 1 + z+/z- + z-/z+ 

2 ^^-'^ 1/2/ , 1/2/ 
^ Z^ / Z- + z_ /z+ 



[St (37) 

(38) 
(39) 



The numerical cancelation of [5]^ fails either when A+z_|_ w 
X-Z- or when A'^X±z± << 1 and these two cases must 
be dealt with separately. In the former case, first order 
Taylor developments of S around the mid value A'^{X+z+ + 
X-Z-)/2 are sufficient whereas in the latter case (typically 
when A'^X±z± « 10"^), series expansions of S in zero 
can be used (Eq. [26l). In both cases, the difference is pro- 
portional to [Az]^ which must be computed as: 



[Az] + 



z+ + z_ 



A+ + A_ 



(40) 



for low photon energy {ujq << 1). With the additions of this 
subsection, Eq. [25] is accurately evaluated for any particle 
and photon energy. 

1.4. Comparisons to other formulae 

This expression was checked extensively against previous 
formulae. The exact differential cross section was first com- 
pared to the original Jones' formula (identical to [T^ when 
the misprints are corrected) and to the formula by NP943 
(Eq. 8.1.7-8.1.22). The latter was derived by integrating 
Eq. [5] in a different order and behaves numerically better 
than the Jones' cross section for Compton up-scattering of 
soft photons (sec Fig. [2|) . Results show perfect agreement 
of the three formulae when the numerical evaluation of the 
three cross sections is accurate. The first moments of the 
distribution: 



cro(po, Wo) 
CTl(po,Wo) 
CT2(P0, Wo) 



da 

——duo 

dio 

[bJ — bjQ)—du) 

duj 



{ll> - UJo} 



I da 
dio 



-dhJ 



(41) 
(42) 
(43) 



- namely the total cross section, the mean energy and the 
dispersion of the scattered photon energy - were also com- 
puted numerically by integrating over the scattered distri- 
bution and they were compared to analytical equation^! 
(NP94, Eq. 3.3.1-3.1.10). Examples are shown in Fig. [3 
Again, good agreement is obtained. Small differences are 
observed but they result from numerical errors in the inte- 
gration of the scattered distribution. 

Compared to previous ones, the cross section derived 
here is found to be accurate over a wider energy range. An 
example of the scattered distribution is shown in Fig. ID As 
summ ed-up in Fig. [2l the original cross section by Jonei 
(fT96l is inaccurate in typical astrophysical situations in- 
volving low energy photons up-scattered by high energy 
particles. The cross section derived by NP94 is accurate in 
most astrophysical cases. When computed directly from the 
published equations, it fails, however, for very high energy 
particles or very high energy photons. Although these cases 
may not be physically relevant since the Klein-Nishina cross 



The equations 
too. In particular; 

I - [X - X^f/Drr, 

Hi = [Ho - AhAi{h. 



in NP94 contain some misprints 
Eq. 8.1.10 should read /i+ — 
/x/xi and Eq. 8.1.22 should read 
)+H/A{h.)A{h+)]/h+ 



A few misprints were corrected. In Eq. 3.3.5, 'I'i4 should read 
(2-i/)_ii — -i/j-io — '!/'-i2)/Ci the last expression for Eq. 3.3.9 should 
be + In^ (2^)/2 - 3.(1/2^), and in Eq. 3.3.10, Voo should 



read 



3/(80 [5(0 - 2/C + (1/2 + 2/e + l/e)k - fi«/2 - 3/2] . 
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Fig. 3. First three moments of the scattered distribution 
as a function of the energy of the incoming photon: the to- 
tal cross section, the mean energy of the scattered photon 
and the dispersion about this mean value. The solid and 
dashed lines show the numerical integration of the expres- 
sion derived in this paper and analytical results given in 
NP94 respectively. The 3 curves are for po = 10~^, 10^ and 
lO**. For each panel, relative errors are also shown. 



section drops at these high energies, numerical issues can 
generate infinite values, propagate and affect numerical so- 
lutions. Nevertheless, it must be noted that a few changes 
in the way some quantities are computed enable good ac- 



curacy over a much larger domain (J. Poutancn, private 
communication) . The cross section derived here is accurate 
for all photon and particle energies and overcomes these 
numerical issues. 



Po 



.Oe + 00 



cj„= l.Oe-05 




10-= 



Fig. 4. The scattered distribution for soft photons {hi/ = 
lO^^nieC^) and mild relativistic particles {p — 1). Solid and 
dotted lines show the distributions obtained with Equation 
[28l and with Jones' formula respectively. 



2. Solving the kinetic equations 

Once the Compton cross section is computed accurately, 
the Compton contribution to the evolution of interacting 
particles and photons is described by the set of equations: 



c^{Po,^o;^^{p)) {po ) dN^ ( Wo ) 
dp 

-Ne{p) / caQ{uJo,p)dN^{ujQ) , (44) 



c^{po,cJo;uj)dNe{po)dN^{ujo) 

-N^iuj) [ cao{Lu,po)dN,{po) , (45) 



dtN,{p) 



where uj{p) = + 7o — 7 is the energy of the scattered 
photon and d(7{po,u!o;i^{p)) = da{po,LL!Q;uj) is the differ- 
ential cross section of the interaction (po,wo) — > (p, w). In 
this section we show that computing these integrals nu- 
merically can also lead to accuracy issues. We present an 
alternative approach that combines this integration with a 
Fokker-Planck treatment and overcomes most numerical is- 
sue s. The basic idea of such a combination was proposed 
by iNavakshin fc Melial (|l998l ) for particles only but was 
not investigated in detail. A more precise analy sis of this 
treatment was presented in iBelmont et al.l ()2008[ ) (see also 
Vurm & Poutanen 200^. Here we investigate the regimes 
where the various methods are valid and we focus on their 
numerical accuracy. 
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2.1. Integral approach 

The simplest way to compute the time evolution of the 
lepton and photon populations is to discretize their distri- 
butions in energy bins. Integration over energy is then per- 
formed bin by bin by summing the contributions from all 
other bins. This process is time consuming since a double 
integral must be performed at each time step. Moreover it 
can be inaccurate in many situations: when the width of the 
scattered distribution is smaller than the bin size, the en- 
ergy change of the scattered particles (or equivalently pho- 
tons) is too small to be resolved by the numerical scheme. 
The particles (or photons) do not scatter far enough in the 
energy space to reach other bins and there is no numerical 
evolution of the distribution. It is well known that most 
numerical issues arise from what happens below the grid 
scale. In this peculiar case however, particles (or photons) 
can undergo many scattering events per unit time when the 
density of scattering centres is large. Errors associated with 
individual scattering events add and they can lead to accu- 
racy issues. This is typically what happens in simulations of 
X-ray binaries: the soft disc photons are up-scattered by the 
high energy particles constituting the corona. This interac- 
tion is known to efficiently cool the corona down. Whereas 
the photons experience large energy changes, the relative 
variation of the particle energy is very small and their cool- 
ing may not be captured by a discrete integral with finite 
resolution. The energy is then not conserved since the pho- 
ton field gains energy while the particle population does 
not lose anjQ. 

The accuracy of this integral method depends on the 
comparison of the energy bin size to the width of the scat- 
tered distribution. The latter depends on the momentum po 
of the incoming particle and the energy wq of the incoming 
photon. Most astrophysical cases involve large ranges of en- 
ergy (typically many orders of magnitudes) which implies 
the need to use logarithmic grids, the resolution of which is 
also energy-dependent. The accuracy of the numerical in- 
tegration thus strongly depends on the region in the lepton 
momentum-photon energy space (poi ^o) the simulation has 
to deal with. In this subsection, the regimes where the inte- 
gral approach remains accurate are investigated. Hereafter 
the subscript is dropped for simplicity. 



2.1.1. Photons 

Numerical integration for the evolution of the photon pop- 
ulation is accurate in the limit where the width of the scat- 
tered distribution Aijj(p, a;) — cjniax(p, w) — Wi„in(p, w) is 
much larger than the local photon bin size Soj. For lin- 
ear grids, requiring that the scattered photon distribu- 



* We note that the energy is not conserved only with non- 
linear grids. For linear grids and identical resolution in the pho- 
ton and electron energy grids {ujo and 70 respectively), the en- 
ergy is balanced to machine precision. In such case indeed, the 
bin size is uniform. As the energy width of the scattered distri- 
bution is by definition the same for both species, if the scattering 
of particles to neighbour energy bins is not captured by the nu- 
merical integration, the scattering of photon is not either, so 
that the integral approach fails simultaneously for both species 
and there is no net energy exchange between the two species. 
Although the integral approach is clearly not accurate in this 
case, the energy balance is computed to machine precision. 



tion spans at least n bins would simply reacH: Aw > 
(n — l)(5w. For logarithmic grids, the same condition reads: 
'^max/'^min > whcrc r^^ = UJ^~^^ / Lj^ is the logarithmic 

increment of the photon grid. By denoting R^j the number 
of bins per decade energy, i.e. the resolution of the pho- 
ton grid, and ei.cj — IQ^^^^^)/^'^ — 1, the condition for the 
integral approach to be valid reads: 



Aw(p, w)/Wniin(p,^) > f-l.a 



(46) 



When the resolution is large enough: ei^^ « 2{n — 1)/Ri^. 
However, as n is typically of several, [n — l)/Ruj can be 
larger than unity for low resolution runs and the exact re- 
lation must be used. The condition is satisfied in a well 
defined region of the (p, w) space. Using Eq. [19] and [20] to 
compute the scattered distribution width, one can write 
an equation in to for the boundary of this region. It is the 
solution of a 2nd order polynomial and the direct numer- 
ical integration is found to be accurate for all photon and 
particle energies except in the region where: 

<p < Pi and min [wf (p), Wj"(p)] < uj < uj^{p) (47) 

with 



Pi 
, ,± 



(£Lc. + 4)-f(ei,^-4)^l + £i, 
4ei,c^ 

i[2(l+p-27) + ei,^(7-p) 



-7 + 



'4(7 -p- l) + eL^(l-2p7-h2p2 
e/2 



(7 - p)e - 2p ■ 



(48) 



(49) 



(50) 



The inverse equation pi(u;) for this boundary would be more 
practical since it would be mono-valued and could be used 
directly in the combined approach (see hereafter). It is the 
solution of a high order polynomial and is most easily found 
numerically. 

Boundaries for the integral approach are shown in Fig. [5] 
for R^ « 9200, 930, 97, 13, 3.8, 2.0, 1.3 and n = 5. For com- 
parison typical runs have a resolution R^ = 1 — 30. In the 
case Ri^ = 13, only the scattering of low energy photons 
{hi' < 200 keV) off sub-relativistic particles {p < 0.2) is in- 
accurately described by the integral approach. Although it 
is not general, the integral method for the photon equation 
can thus be used to address many of astrophysical situa- 
tions, such as those involving high energy particles. 

2.1.2. Particles 

Identically, the integral approach for particles is valid 
when the width of the scattered distribution is larger 
than the particle bin size. For logarithmic grids this 



„n-l 



where p 



max /min 



condition is: Pmax/Pmin > 

•^(7 -t- w — Wijiin/max)^ ~ 1 are the largest and smallest mo- 
menta of the scattered particle distribution and Vp = 
p^+^/p* is the logarithmic increment of the particle grid. 
By defining Rp and ei.p as for the photons, it yields: 

Ap(p,w)/Pmin(p,w) > ei,p . (51) 



^ The —1 comes from the fact that when the width of the 
scattered distribution is exactly nSui, it is discretized in n + 1 
bins in all situations except when it is exactly centred at the 
centre of a bin. 
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Fig. 5. Photon equation: boundaries V\{'^) for the 
integral approach (solid lines) and pfpC"^) for the 
Fokker-Planck approximation (dashed lines) for 
ei,„ = 10-3, 10-2, 10^1, 1,10, 102,103 and epp = 
10-3,10-2,10-^,1 respectively. The integral approach 
and the FP approximation are valid respectively right and 
left of the related boundaries 
(see Eq.ISI]). 



The dotted line shows oJq (p) 




Fig. 6. Particle equation: boundaries aji(p) for the 
integral approach (solid lines) and uj-p^(j)) for the 
Fokker-Planck approximation (dashed lines) for 
ei^p = 10-3,10-2,10-1,1,10,10^,103 and epp = 
10-3,10-2,10-1,1 respectively. The integral approach 
and the FP approximation are valid respectively above 
and below the related boundaries. The dotted line shows 



Contrary to the photon case, there is no simple equation 
for the boundary defining this region and the implicit equa- 
tion [51] must be solved numerically. The boundary Ct'i(p) is 
computed easily since there is only one positive solution to 
the equation on lo and it happens to be always smaller than 
Wq(p). a simple approximation for this boundary, accurate 
in all regimes to better than 30%, is: 

c.i(p) « . (52) 

2p2 + (1 + ej p - ^1 + ei,p)p -f ei^p -t- 2 

When the resolution is good enough [R^ >> 1, ei,p << 
1), condition [SD reduc es to that used in Eq. 78 of 
IVurm fc PoutanenI (j2009( ) for electron down-scattering: > 
ei,p/(47) for relativistic particles. 

Fig.[6]shows this boundary for various resolutions {Rp w 
9200,930,97,13,3.8,2.0,1.3 and n = 5). Contrary to the 
equation for photons, solving the evolution equation for 
particles with a simple numerical integration requires very 
high resolution to guarantee a correct accuracy. For exam- 
ple, soft photons of energy ui = 10~^ scattering off high 
energy particles (7 = 102) are accurately described by the 
integral approach only if Rp > 900, which is far beyond 
current desktop computer capabilities and the constraint 
is even more severe when mildly relativistic particles are 
involved. 

2.1.3. Accuracy 

If conditions [46] and [51] are not satisfied over the entire 
simulation domain, numerical computation of the Compton 
scattering may lead to inaccuracy. In particular, energy is 
not conserved when logarithmic grids are used (as explained 
before), and its measure provides a good way to estimate 
numerical error^. Fig. [7] shows the error on energy when 

^ For linear grids, other accuracy indicators should be used. 



soft mono-energetic photons are up-scattered by high en- 
ergy particles. The numeri cal computation was performed 
with the code presented in iBelmont et al.l (|2008) , on a sin- 
gle time step, using an explicit scheme on logarithmic en- 
ergy grids, and turning off all processes/injection/escape 
but Compton scattering. The contribution of Compton 
scattering was forced to be computed by the integral ap- 
proach for the equation on particles. The contribution to 
the equation on photons was also computed with the in- 
tegral approach but given the energy range and the grid 
resolution considered here, this method is accurate in this 
case. Both initial distributions were set to zero except in 
one energy bin, and they were normalised to unity to keep 
the number of scattering events constant as the photon en- 
ergy is varied. The error was measured by computing the 
total energy lost by particles dtEp and comparing it to the 
energy gained by photons dtE^. 

For low energy photons, the scattering takes place far 
below the boundary in the (p, uj) plane of Fig. [S] The scat- 
tered distribution typically spans less than 2 bins. The en- 
ergy change of particles is not captured at all by the numer- 
ical scheme and the error is 100%. As the photon energy 
increases, the scattered distribution widens and the error 
decreases. From the cases presented here, we find that an 
error less than 1% corresponds to scattered distributions 
that span more than 5 bins. When n > 10, the error is 
dominated by other weak numerical errors. 

By the choice of particle and photon energies, we have 
only focused here on the error made in the equation for 
particles. A similar study can be done for photons. However, 
it would lead to very similar results and contrary to the case 
of particles, integration errors for photons are less relevant 
to astrophysical applications. 
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Fig. 7. Numerical error on energy conservation, for a sin- 
gle time step, as a function of the incident photon energy 
ujQ, for particles of momentum po = 10^, 10'^, 10"*. In the 
case pq — 10'^, the boundaries between zones where the 
scattered distribution spans 2,3,4,5,6,7,8,9, and at least 10 
energy bins are over-plotted (vertical dashed lines). The 
error is computed as the relative error between the total 
energy losses of the photon population and the total en- 
ergy losses of the electron population: SE/E — {\dtEp\ — 
\dtE^\) / {\dtEp \ + \dtE^\). For this run the resolutions are: 
Rp = 64/6 and R^, = 256/12. 



2.2. Fokker-Planck approach 

Alternatively, Compton scattering can be described by the 
Fokker-Planck method. This approximation assumes that 
the cross section and the initial particle (or photon) dis- 
tribution are only weakly dependent on the energy of the 
incoming particle (or photon) on the scale of the scattered 
distribution width. When the width of the scattered dis- 
tribution is small, a second order Taylor expansion of the 
integral equation can be performed, which gives the well- 
known Fokker-Planck equations: 



dtNe 



2 \P \P 



with 

Auj , — 



4 • n 



cCTi;2(p; w)dA^e(p) 



(53) 
(54) 

(55) 
(56) 



As for the integral approach, the validity of the FP ap- 
proximation depends on the energy of the incoming particle 
and photon. Although the real criteria should in principle 
depend on the initial distributions and should be computed 
at each time step of a time dependant simulation, it is more 
convenient to define an approximate but general criteria. 



width A7 of the scattered distribution. There is no unique 
way to define such a condition. By assuming that the typical 
energy scale for the variation of this quantity is the initial 
kinetic energy: 70 — 1, then the FP approach is found to 
give a correct description of the particle evolution if: 



A7(w,7)/(7 - 1) < epp 



(57) 



where epp << 1 is a small parameter that describes the 
accuracy of the FP approximation: the smaller epp the more 
valid the FP approach. The equation for the boundary of 
this region of the (p, to) space reads: 



U! < CJpp(p) 

where 




- epp(l - 1/7)/? 
epp(l - l/7)//3 



1 



(58) 



(59) 



In the limit epp << 1, a simple approximation for all par- 
ticle energy is: 



UJpp 



epp/? 
4(7 + 1) 



(60) 



which reduces to LOpp^p ^ epp/ (47) in the ultra-relativistic 
regime (p » 1). Boundaries cjpp(p) are shown in Fig. [S] for 
epp = 10-3, 10-^ 10'^ and 1. The Fokker-Planck approxi- 
mation for particles is a good approximation when very low 
energy photons are up-scattered by mid-relativistic parti- 
cles. It fails easily when the particle energy becomes large. 
For example, by setting epp =0.1, the scattering of pho- 
tons off particles with energy p = 10'^ can only be described 
by the FP method when the photon energy is less than 
uj = 10"^. However, when one does not need such a preci- 
sion and sets epp ^ 1, the FP method is valid over a large 
range of energy. 

2.2.2. Photons 

Similarly, the FP method for photons assumes that the 
quantity N^{uJo)d(7{'yo,'^Oi^) does not vary much with uiq 
in the width Auj of the scattered distribution. By assuming 
that the typical energy scale for the variation of this quan- 
tity is of the order of the incoming photon energy loq, the 
FP approach for photons is valid if: 



Auj{u!, 'y)/uj < epp 



(61) 



When solved as a function of lo, the equation for the bound- 
ary is again the solution of a second order polynomial and 
the condition [HI] is satisfied in the region defined by: 



P <Pfp 
where 

Pfp 

UJpp 



and min [ti)pp(p), Wpp (p)] <uj< uj+p{p) (62) 



(4-epp)-(4-f epp)Vr 



epp 



2(1 -p 



4epp 

7) + epp (7 - 



P) 



7 - p) 4- epp(H- 2p7 2p2 



(63) 



(64) 



2.2.1. Particles 

The FP method for particles assumes that the quantity 
Ne{'yo)da{"fQ,ujQ;uj) does not vary much with 70 in the 



Wpp 



epp 



cpp 



/Vp' 



1 + epp//3 + V(l + epp/Py - e^p/p2 



(65) 
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As for the integral approach, the inverse equation for this 
boundary pfp(lu) is more practical to use and is best eval- 
uated numerically. 

Boundaries ppp are plotted in Fig. [5] for epp — 
10^^, 10^^, 10"-"^, and 1. Contrary to the equation for par- 
ticles, the photon evolution is only poorly described by the 
FP approximation and only low energy photons (w < 0.1) 
scattering off low energy particles {p < 0.1) undergo small 
angle scattering that can be accurately described by the FP 
approximation. 

2.3. Combined approach 

As can be seen in Fig. [5] and [6l the two methods happen 
luckily to give a correct description of Compton scattering 
in regions of the {p, uj) space that can be complementary. 
This suggests that combining the two methods is a good 
way to overcome numerical accuracy issues in computing 
the effects of Compton scattering. However, this can only 
be done if the validity regions for the two methods cover 
the entire simulation domain, that is, for given ei^p, ei^cj and 
epp, if: 

ppp(a)) > pi(uj) , (66) 

i^Fp{p) > i^i{p) Vp (67) 

for the equation on photons and particles respectively. 
Typically, this implies that ei^^i < epp and ei^p < epp, which 
sets a constraint on the grid resolution: 

1 



epp 



(68) 



As good accuracy for each methods requires n > 5 and 
epp << 1, this puts strong constraints on the resolution. 
For example, by setting epp — 0.1 and n = 5 this leads to 
i? > 40. For particle grids that typically span 5 orders of 
magnitude, it implies the need to use 200-bin grids accessi- 
ble to most desktop computers. However, for photon grids 
that typically span 15 orders of magnitude, it corresponds 
to 600-bin grids, which requires high computing power. 

For each equation (for photons and for particles), one 
can define an average boundary when the validity regions 
of both methods overlap: 



LOcip) = (cjicjpp)^^^ and pdu)) = {piPfp)^^^ 



(69) 



Finally, a combination of the two methods is achieved by 
solving the following equations: 



dtN,{p) = 



dp ^A^<^-N, 



2 " Vp '^Vp 

/•CO ^- 

dNe{po) / c—{po,ojQ;uj{p))dN^{ujQ) 
^c(po) dp 



- Ne{p) / cao{too,p)dN^{ujo) , (70) 
dtNUco) = {Al<P^N^) + \dl. {Dl<^^N^) (71) 
dN^{LJo) / c— (po,^o; w)dA^e(po) 



N^{ui) / cao{uj,po)dNe{po) 



(72) 



where integrations in the integral part are performed only 
above the boundary energy lOc and momentum pc and where 
the FP coefficients are computed as integrals of the cross 
section only below these boundaries: 



4 ■ n — 



Pc(w) 



ccri.^2ip;uj)dNe{p) 



■'c(p) 



(73) 
(74) 



Ap;Dp = / cai.^2ip;^)dN^{uj) . 
Jo 

Fig. [5] shows the steady photon and particle distributions 
for a more realistic simulation. Particles are injected with 
a mono-energetic distribution at high energy (7 = 10^) 
at a constant rate (the compactness of which is set to 
^in,c = (TTP/imeC^R) = 10, where P is the injected power) 
and escape at a constant rate with the speed of light to al- 
low for a steady state. Soft photons are also injected with a 
black body spectrum of temperature ksT = 3 x 10~^TOeC^ 
and a fl ux of compactness / jn.^ — ut L soft / {rrieC^R) — 1 
(see e.g. 'Bel mont et al.ll2b08l for detailed definitions of the 
compactness parameters). We start with an empty sys- 
tem and the code simultaneously evolves both distributions 
with time until they reach a steady state. The only process 
turned on is Compton scattering. Fig. [5] shows the results 
of three kinds of runs where the contribution of Compton 
scattering to the evolution of the particle distribution was 
forced to be treated with the integral (several resolutions), 
the Fokker-Planck and the combined approaches. The dis- 
tributions look very different. 

As expected, the integral approach fails to capture effi- 
ciently the particle cooling for low resolution runs: the lower 
the resolution, the steeper the slope of the high energy tail 
of the particle distribution. These excess high energy par- 
ticles more efficiently up-scatter the soft photons and the 
resulting spectra are harder. In steady state, the total en- 
ergy loss must balance the total injected power. Power is 
supplied through particles h-n.e = 10 and through photons: 
^in = ^in,G + Un,ui = H- Energy is also lost both through 
particles and photons. The steady particle distribution has 
low density t = ctR J dN^ — 2.3 x 10~^ so that the energy 
lost through particle escapes is rather small (^out.o ~ 1) 
and most of the energy escape as radiation. As numerical 
errors prevent particles from cooling as they should, the 
overall system actually gains energy artificially and the to- 
tal energy loses are measured to be larger than the injected 



power: Zout = I 



out,aj ~t- ^out,e 



19.6, 14.2, 12.2, and 11.4 from 



lower to higher grid resolution. The corresponding effec- 
tive electron temperatu res are fcgTe = 6.23, 4.68, 3.95, and 
3.60 ruec^ (see Eq. 2.8, lCoDpilfT99l . 

The Fokker-Planck equation provides good energy con- 
servation and the plasma luminosity balances the injected 
power: ^out = 11.0. However, energy conservation does not 
guarantee accurate results and significant deviation is ob- 
served in the low energy part of the particle distribution, 
where the Fokker-Planck approximation is expected to fail. 
In the particular case presented here, this has no effect on 
the emitted spectrum since the spectrum is dominated by 
the Comptonisation by high energy particles. Nevertheless 
it could strongly affect the predicted spectrum when other 
processes such as Coulomb collisions or pair annihilation 
are involved. The electron temperature is measured to be: 
/c^Te = 3.40 rnec^. 

The combined approach produces the best results. With 
low grid resolution {Rp = 64/6), it gives the correct photon 
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Fig. 8. Particle distribution and spectra in steady state 
when Compton scattering for particles is treated with 
the integral method for different resolutions: Rp = 
512/6,256/6,128/6, and 64/6 (solid lines), when it is 
treated with the combined approach with epp = 0.1 and 
n = 10 (dashed line), and with the Fokker-Planck approx- 
imation (dotted line). For these runs, the power injected 
in particles and photons is: ^in,c = 10, ^m.u = 1 and the 
source size is set to i? = 10^ cm. The temperature of the 



seed photons is ksT /me 



3 X 10 and the energy of 



the injected particles is 7 = 10^. For the combined and FP 
approaches, the particle resolution is: Rp = 64/6. For all 
runs, the resolution of the photon grid is R^ = 512/9. 



spectrum and the correct particle distribution, it conserves 
energy to better than 1% accuracy (Zout = 11.0) and it 
was checked that the shape of the steady distributions does 
not depend on the grid resolution (not shown here). The 
electron temperature is fc^Te = 3.39 mgC^. In comparison, 
the integral approach requires more than 10 times higher 
resolution to produce similar results, which also means at 
least 10 times more computational time. Extensive tests 
show that the combined approach is very efficient in most 
cases. It takes advantage of the two approaches. The way 
the two methods are associated is numerically simple to im- 
plement and it can easily be shown that the total number 
of particles and photons as well as the the total energy are 
conserved. Moreover even when conditions [66] and [62] are 
not exactly satisfied this method still provides good accu- 



racy in many cases since, as long as these conditions are not 
strongly violated, only a small faction of particles and/or 
photons is treated inaccurately. Also, even if the evolution 
of a large fraction of particles and/or photons is not de- 
scribed correctly by the numerical scheme they might have 
only a minor contribution to the overall distribution evo- 
lution in some specific cases. For these reasons some simu- 
lations might reach a correct final accuracy, despite violat- 
ing the accuracy condition. However, this is very problem- 
dependent and the accuracy should be checked very care- 
fully when these conditions are not satisfied. 

Conclusion 

In this paper we have investigated numerical issues related 
to the computation of isotropic Compton scattering in nu- 
merical codes. We have derived a form of the distribution 
of Compton scattered photons that is free of numerical can- 
cellations. This form is exact with no limitation on the pho- 
ton and electron energies. The numerical accuracy resulting 
from the direct use of the cross section has been studied for 
various grid resolutions and equations have been proposed 
for the boundaries where this method must be replaced 
by an approximate Fokker-Planck treatment to guarantee 
sufficient accuracy. All results presented here have been in- 
cluded i n and extensively test ed with the new code devel- 
oped bv lBelmont et al.l ((2OO8O . 
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